Half-quantum vortex state in a spin-orbit coupled Bose-Einstein condensate 
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We investigate theoretically the condensate state and collective excitations of a two-component 
Bose gas in two-dimensional harmonic traps subject to isotropic Rashba spin-orbit coupling. In the 
weakly interacting regime when the inter-species interaction is larger than the intra-species interac- 
tion (g^i > g), we find that the condensate ground state has a half-quantum-angular-momentum vor- 
tex configuration with spatial rotational symmetry and skyrmion-type spin texture. Upon increasing 
the interatomic interaction beyond a threshold g c , the ground state starts to involve higher-order 
angular momentum components and thus breaks the rotational symmetry. In the case of g-fi < g, 
the condensate becomes unstable towards the superposition of two degenerate half-quantum vortex 
states. Both instabilities (at g > g c and g^i < g) can be determined by solving the Bogoliubov 
equations for collective density oscillations of the half-quantum vortex state, and by analyzing the 
softening of mode frequencies. We present the phase diagram as functions of the interatomic inter- 
actions and the spin-orbit coupling. In addition, we directly simulate the time-dependent Gross- 
Pitaevskii equation to examine the dynamical properties of the system. Finally, we investigate the 
stability of the half-quantum vortex state against both the trap anisotropy and anisotropy in the 
spin-orbit coupling term. 

PACS numbers: 05.30.Jp, 03.75.Mn, 67.85.Fg, 67.85.Jk 



I. INTRODUCTION 

Owing to the unprecedented control in interatomic in- 
teraction, geometry and purity, atomic quantum gases 
have proven to be an ideal many-body platform for 
exploring fundamental quantum states, such as Bose- 
Einstein condensates (BEC) [I], strongly interacting uni- 
tary Fermi superfluids O [3] and Mott-insulating states 
@]. One of the latest achievement concerns the spin- 
orbit (SO) coupling in an ultracold spinor Bose gas of 
87 Rb atoms [5] , induced by the so-called "synthetic non- 
Abelian gauge fields". Novel quantum states may be an- 
ticipated in the presence of SO coupling [fJHTB], Indeed, 
for a homogeneous SO coupled spin- 1/2 Bose gas with 
intra- and inter-species interactions (g and g fl ), a sin- 
gle plane-wave or a density-stripe condensate state has 
been predicted [8\, depending on whether g is smaller or 
larger than <j> . Interesting density patterns have been 
observed in the theoretical simulations for an SO cou- 
pled spinor condensate, in the absence jSJ [TTJ []2J [T71 [TS] 
or presence |14fJTf5] of rotation. The phenomenon of self- 
trapped BECs has also been proposed, in particular, in 
one-dimensional (ID) geometry [7]. 

In this work, we show that in a Rashba SO coupled, 
weakly interacting spin- 1/2 Bose gas in two-dimensional 
(2D) harmonic traps, all bosons may condense into a 
non-trivial half-integer angular momentum state (or a 
half-quantum vortex state) with a skyrmion-type spin 
texture. We solve the mean- field Gross-Pitaevskii equa- 
tion (GPE) for its density distributions and spin textures, 
and obtain its collective excitation spectrum by solving 
the Bogoliubov equation and by directly simulating real- 
time propagation of the GPE ground state under per- 
turbation. The condensation of an SO coupled spin-1/2 



Bose gas into a half-quantum vortex configuration was 
first suggested by Congjun Wu and co-workers in 2008 
and its existence was discussed under the condition that 
the interaction is SU(2) symmetric, i.e., g = g fl |13| . 
Here, we explore systematically the parameter space for 
the half-quantum vortex state and analyze its stability. 
We present a phase diagram for the half-quantum vortex 
state as functions of the SO coupling and the interatomic 
interaction strengths. We also investigate the dynamical 
properties of the half-quantum vortex state by directly 
simulating the time-dependent GPE. Finally, the stabil- 
ity of the half-quantum vortex state against both the trap 
anisotropy and anisotropy in the spin-orbit coupling term 
is examined. 

Our main results are summarized in Fig. [T] The 
half-quantum vortex state (the phase I) is the ground 
state if the intra-species interaction is smaller than the 
inter-species interaction ( g < g., ) and if the interaction 
strength is below a threshold (g < g c ). Otherwise, it 
becomes energetically unstable towards a superposition 
state of two degenerate half-quantum vortex states (the 
phase IIA) , or a state involving higher-order angular mo- 
mentum components (the phase IIB). With decreasing 
the dimensionless SO coupling strength Xso, the thresh- 
old g o becomes exponentially large, leading to a large 
parameter space for the half-quantum vortex state (see 
Fig. 14 1. It is therefore feasible to be observed in the cur- 



rent experiments with ultracold SO coupled spinor Bose 
gases of 87 Rb atoms. 

The rest of the paper is organized as follows. In the 
next section, we outline the model Hamiltonian and dis- 
cuss briefly the existence of half-quantum vortex state 
in the non-interacting limit. In Sec. Ill, we present the 
numerical procedure of solving the GPE and Bogoliubov 
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Figure 1: (color online). Phase diagram at two dimension- 
less SO coupling strengths, Xso = 1 (a) and Xso = 4 (b). 
The half-quantum vortex state (the phase I) becomes un- 
stable when the intra-species interaction is larger than the 
inter-species interaction (g > g n , the phase IIA) or when 
the interatomic interactions are sufficient strong (g > g c , the 
phase IIB). The insets shows the density patterns of the spin- 
up and spin-down bosons in the phases I and IIA. We note 
that, the critical interaction strength g c increases rapidly with 
decreasing the SO coupling strength Xso- 



equations and discuss the typical density distributions 
and collective mode behaviors of the half-quantum vortex 
state. The collective excitation spectrum obtained from 
the Bogoliubov equation is compared to a direct simula- 
tion of the time-dependent GPE. In Sec. IV, we analyze 
the stability of the half-quantum vortex state by mon- 
itoring the softening of collective mode frequencies and 
by comparing the energy with that of some competing 
states. The phase diagram is then constructed as func- 
tions of interatomic interactions and SO coupling. The 
stability against the anisotropy in trapping potential and 
in spin-orbit coupling term is also carefully examined. Fi- 
nally, we summarize in Sec. V and give some concluding 



remarks. 



II. THEORETICAL FRAMEWORK 

We consider a two-component Bose gas confined in 
a 2D isotropic harmonic trap V(p) = Mlo\(x 2 + 
y 2 )/2 = Muj]_p 2 /2 with a Rashba SO coupling Vso — 
—iXni^xdy — a y d x ), where is the Rashba SO cou- 
pling strength and a x , a y , and a z are the 2x2 Pauli 
matrices. The model Hamiltonian H — J dr[Ho + Hi n t] is 
given by, 
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where r = (x,y) and vj/ = [vj/^.(r), ^^(r)] T denotes the 
spinor Bose field operators in a collective way, and the 
chemical potential \i is to be determined by the total 
number of bosons iV", i.e., J dr^f^ = N. For simplicity, 
we have assumed equal intra-species interaction strength 
fftt = 9ii = 9- I n experiments, the two-dimensionality 
can be readily realized by imposing a strong harmonic 
potential V(z) — Mu> 2 z 2 /2 along axial direction, in 
such a way that p,ksT <C huj z [TH]. For the realis- 
tic case of 87 Rb atoms, the interaction strengths can 
be calculated from the two s-wave scattering lengths 
a ~ lOOas and Ofj., using g = v / 87r(/i 2 /M)(a/a z ) and 

9U = 



87r(fi 2 /M)(a^i/a z ), respectively. Here 



\Jhj (A1uj z ) is the characteristic oscillator length in z- 
direction. 

For a weakly interacting Bose gas at zero tempera- 
ture, we assume that all the bosons condense into a 
single quantum state ^(r) = [4>-|-(r), $^(r)] T . Following 
the standard mean-field theory [5U], we separate the 
field operator into a condensate and a fluctuation part, 
^(r) =$ cr (r)+^ cr (r). Keeping up to the quadratic 
terms in ^o-(r), this separation leads to H = J dr[Hop + 
Hr], where the condensate part is given by, 
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and the fluctuation part Ht = ^Hsog 1 ^ with 



H-Bos — 



H st +g\^ V so + g n ^l 

Kl + .9tl^l H Sl +g\^\ 2 

9 ($}) 9n^l 

9U®t®t 9 



H 



9^ 

gn®t®i 

9\®t 



9*1 

-Vl + gn^i 



(4) 



V so + g n ^l Hs.+g^xY 



3 



Here H osc = -H 2 V 2 /(2M)+V (p), U s ^ = H osc +g |$ t l + 

g n |$J 2 - ii and U H = H osc + g n |$ t | 2 + g |*J 2 - /x, 

V so = -i\ R (d y + id x ) and = -iX R (d y - id x ), 

and we have introduced a 4 x 4 Nambu spinor — 
[* t (r),^(r),^(r),^ (r) ]T. ^ 

The condensate wave-function can be obtained from 



the GP equations 6HGp/8®{r) 



-i\ R (d y - id x ) 



-iX R {d y 

flsi 



= o EH 

id x ) 



or explicitly. 



*t ( r ) 
(r) 



= 0. 



(5) 

At zero temperature, we assume a single condensate state 
with zero quantum depletion, so that the condensate 
wave- function is normalized by Jc?r[|$-|-| 2 + |$^| 2 ] = N, 
where N is the total number of bosons. The equation 
becomes simplified if we write <I>^ = N 1 / 2 ^ and <£>^ = 
N 1 ' 2 ^ and use accordingly the interaction strengths 
g(N — 1) and g^(N — 1). The normalization condition 

becomes J dr [\(j>^\ 2 + |0^| 2 ] = 1. 

The quasi-particle wave-functions with energy Hu sat- 
isfy the Bogoliubov equations |20| . 
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and is normalized by Jo?r[|u-j-| +|itj,| — |uf| — |i>4.| ] = 1. 
These Bogoliubov quasi-particles correspond to the dif- 
ferent collective density oscillation modes around the con- 
densate with the frequency u [5T]. It is easy to see that 
the wave-function [v* (r) , uf (r) , ut (r) , u| (r)] T is also a 
solution of Eq. ([6]), but with energy — Hu. This is antic- 
ipated for the usual Bogoliubov transformation. Phys- 
ically, we should restrict to a non-negative mode fre- 
quency, u > 0. 

In harmonic traps, it is natural to use the trap units, 
i.e. to take Hu± as the unit for energy and the harmonic 
oscillator length aj_ = WH/ (Mu±) as the unit for length. 
This is equivalent to set H = ks = M = u± = 1. For the 
SO coupling, we introduce an SO coupling length a\ — 
H 2 /(MX R ) and consequently define a dimensionless SO 
coupling strength Xso = aj_/a\ = (M/H 3 )X R /^/u±. 
In an SO coupled spin-1/2 BEC of 87 Rb atoms as realized 
recently by the NIST group [5:, Xso is about 10. In 
the typical experiment for 2D spin-1/2 87 Rb BECs [15] , 
the interatomic interaction strengths are about g(N — 
1) ~ 9U( N ~ !) = 102 ~ l0 3 Hu J _/a 2 L . These coupling 
strengths, however, can be precisely tuned by properly 
choosing the parameters of the laser fields that lead to 
the harmonic confinement and the SO coupling. 



A. Single-particle solutions 

The appearance of the half-quantum vortex state may 
be easily understood in the non-interacting limit |13| . 



In the absence of interatomic interactions, the single- 
particle wave-function [</)^ (r) , <f>^ (r)] T with energy e is 
given by, 



-iX R (d y - id x ) 



-iX R (d y 



id x ) 



(7) 

In polar coordinates (p,tp), we have —i(d y ± id x ) — 
e Tlip [±d/dp — (i/p)d/d(p]. Because of the isotropic har- 
monic potential V (p), the single-particle wave- function 
may have a well-defined azimuthal angular momentum 
l z = r7i and may take the form, 



<^m(r) 



<h{p) 



/2tt 



(8) 



This state also has a well-defined total angular momen- 
tum j z = l z +s z = 777+1/2. In general, we may denote the 
energy spectrum as e nm , where n = (0, 1, 2...) is the quan- 
tum number for the transverse (radial) direction. There 
is an interesting two-fold degeneracy of the energy spec- 
trum: any eigenstate <j){r) — [0f(r), </^(r)] T is degener- 
ate with its time-reversal partner T</>(r) = (ia y C)(f)(r) = 
[<p^(r), — (/>|(r)] T . Here C is the complex conjugate oper- 
ation. This Kramer doublet is the direct consequence of 
the time-reversal symmetry satisfied by the model Hamil- 
tonian. It preserves as well in the presence of interatomic 
interactions. As a result, we may restrict the quantum 
number m to be non-negative integers, as a negative m 
can always be regarded as the time-reversal partner for 
a state with m > 0. 

To solve numerically the single-particle spectrum, we 
adopt a basis-expansion method. To this end, we expand 
first, 



<h(p) = ^A k R km (p) , 
fc 

hiP) = ^2 B kRkm+l (P) , 



(9) 
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where 
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v |m| 2 

a ± \ (k + \m\)l \a ± J k v a 2 / v ; 



is the radial wave-function of a 2D harmonic oscillator 
Hosc with energy (2k + \m\ + l)Hu±, and is the asso- 
ciated Legendre polynomial. Then, we have the following 
secular matrix, 



n osct M 1 
m n oscl 













B k 


= e 


Bk 



(12) 



where the matrix elements are given by (for m > 0) 

^osct,kk' = Hu ± [2k + m + l}S kk ', 

= Hu± [2k + (m + 1) + 1] 8kk', 

Mkk' = hu±Xso Vk' + m+ 15 kk > + vWSkk'-i 
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Diagonalization of the secular matrix Eq. (12 1 leads 



to the single-particle spectrum and single-particle wave- 
functions. In numerical calculations, it is necessary to 
impose a cut-off fc max for the radial quantum number k 
of the 2D harmonic oscillator. For Xso < 20, we find that 
&max = 256 is already sufficiently large to have an accu- 
rate energy spectrum. With this cut-off, the dimension 
of the secular matrix in Eq. ( 12 1 is 2/c max = 512. 
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Figure 2: (color online), (a) Single-particle energy spectrum 
at Xso = 1- (b) The density profiles for the single-particle 
state with m = at Ago = 1- ( c ) The IF-function for the m = 
single-particle state as a function of SO coupling strength. 
It is always positive at arbitrary SO coupling strength. 

In Fig. 2a, we show the single-particle energy spectrum 
at Xgo = 1- For arbitrary SO interaction strength, we 
find numerically that the doublet single-particle ground 
state always occurs at m = (or m = — 1 for its time- 
reversal partner state). 



Here a and are two arbitrary complex numbers satis- 
fying H 2 + \(3\ 2 = 1. 

In the presence of very weak interatomic interactions 
such that g(N - l)a 2 ± ,g n (N - l)a]_ < Ae, where Ae is 
the energy difference between the single-particle ground 
state 4>q (r) and the first excited state 0i(r), we may de- 
termine the superposition coefficients a and j3 by min- 
imizing the GP energy, £ GP [0 s (r)] = / dr Hgp[<M f )]- 
After a simple algebra, we find that, 



AE = E GF [ct> s (r)}-E GP [<f> (r)], 

= (g n -g)(N-l)\aP\ 2 W[M*)], 
where the IF-function is given by, 



2x2 



W[<f,(r)]= / A-Kl^rl -1^1 ) 



(14) 
(15) 



(16) 



Therefore, a half-quantum vortex state is preferable if 
(.9ti ~ 9)W > 0. Otherwise, an equal-weight superposi- 
tion of two degenerate half-quantum vortex states with 
M = \fi I = Vv2 w iH be the ground state. As shown in 
Fig. 2c, the IF-function for 0o( r ) is positive for arbitrary 
SO coupling. We thus conclude that a half-quantum vor- 
tex state should appear at weak interatomic interactions 
provided that the inter-species interaction is larger than 
the intra-species interaction (g-\\. > g). 



III. DENSITY DISTRIBUTIONS AND 
COLLECTIVE EXCITATIONS 

Let us now consider finite interatomic interactions, 
by solving the GPE for density distributions and spin- 
textures, and the Bogoliubov equation for the collective 
density excitations. 



B. Appearance of the half-quantum vortex state 

The single-particle state with m = 0, <fio(r) = 
[cf>^(p),cf>l(p)e lv ] T /\/2tt, has a half-quantum vortex con- 
figuration [131 I22| . as the spin-up component stays in 
the s-state while the spin-down component in the p- 
state and the resulting spin texture is of skyrmion type 
(see Fig. 2b for density distributions and Sec. IIIB for 
more discussions on spin-texture). In the absence of in- 
teractions, however, there is a degenerate time-reversal 
state, T<j)o(r) = [4>i(p)e~ lip , —cf)^ {p)\ T /y/2TT, which is also 
a half-quantum vortex state. Therefore, in general, the 
ground single-particle state is a superposition of two de- 
generate half-quantum vortex states of (f> (r) and T<t>o ( r ) ■ 
which takes the form <j> s (r) = a(j>o (r) + f3T<fio (r) , or ex- 
plicitly, 



A. GPE solutions of the half-quantum vortex state 

For the half-quantum vortex condensate state with 
m = 0, the GP equation becomes Cgp{4>1; (p) j <t>X (p)] = 
0, where 
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s,0 



b t + 9n<f>l 



x R (d p + l/p) 
(,r " X R (-d p ) H.,i + g n <% + g<% 

(17) 

g = g(N - 1)/(2tt) and g n = g n (N - 1)/(2tt), and 

n s . m = ~[h 2 /(2M)}[d 2 /dp 2 +(i/p)d p -m 2 /p 2 }+v(p)-p. 

The numerical procedure for solving GPE is very similar 
to that for single-particle states in Eq. <\12\ . We expand 
<h(p) = J2k A kRka(p) and <pi(p) = J2k B k R ki (p), and 
obtain the secular matrix (with m — 0), 
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where 



L \,kk' 



X 



pdpRko (p) {g<f% + gn<t>f) R k'o (p) ,( 19 ) 

pdpR kl (p) (g n $ + g<f>l) R V1 (p) .(20) 



The chemical potential is given by the lowest eigenvalue 
of the secular matrix. Due to the non-linear terms of 
Zf,fcfc' and T^^k' > we have to update the condensate wave- 
functions and densities iteratively. To overcome the large 
non-linearity, we use a simple mixing scheme by setting a 
small parameter < 7 < 1 and replace the previous den- 
sity <t>l,old h y ( X _ 7) + where <j>l is the density 
calculated in the current step |23| . The choice of 7 de- 
pends on the interaction strengths. It becomes smaller 
for larger g and g^. We run the iteration until conver- 
gence is achieved within a set tolerance. We have checked 
that this procedure of solving GPE is stable for interac- 
tion strengths up to g(N — l),g^i(N — 1) < 10 3 fiu;^/ 'a 2 ^. 
For even larger non-linearity, it seems to be impractical 
to expand the condensate wave-function using the 2D 
harmonic oscillator basis. Therefore for large interac- 
tion strengths, we use a time-splitting spectral method 
(TSSP) technique to solve the coupled GP equations and 
obtain the ground state by imaginary-time propagation 
[2~31l2"5] . For small interaction strengths, results obtained 
from TSSP are identical to those obtained from the basis- 
expansion method. 



that the density distributions are flattened significantly 
by interatomic interactions, as anticipated. The 2D con- 
tour plot of the spin-up and spin-down density patterns 
of the half-quantum vortex state is shown in the inset of 
Fig. [I] (in the phase I). 

To gain more insights of the half-quantum vortex state, 
it is useful to calculate the spin vector 



V ' 2 |$| 2 



and the skyrmion density 



^skyrmion ( r ) — ~j S • [d x S X cLS] 
47T 



(21) 



(22) 



The skyrmion density is a measure of the winding of the 
spin profile. If it integrates to 1 or —1, a topological 
stable knot exists in the spin texture. 




Figure 4: (color online). Contour plots of the three compo- 
nents of spin vector S (r) at Aso = 1, g(N — 1) = 40fkj±/a 2 L 
and g u /g = 1.1. 



B. Density distributions and spin textures 
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Figure 5: (color online) (a) Two-dimensional vector plot 
of the transverse spin vector (S Xy S y ) at Aso = 1, g{N — 
1) = 40&J x / 'a± and g^i/g = 1.1. The color and length 
of arrows give respectively the orientation and the magni- 
tude of (S Xl S y ). (b) The corresponding skyrmion density 

^skyrmion (l*)- 



Figure 3: (color online). Density distributions at Ago = 1 and 
g(N - 1) = 40/iajx/ai (a) and at Aso = 4 and g(N - 1) = 
huj±/a 2 L (b). Here, the ratio gu/g = 1.1. 

In Fig. [3] we present the radial density distributions of 
the half-quantum vortex condensate state at two SO cou- 
pling strengths: Aso = 1 and Aso = 4. The increase of 
the SO coupling leads to more oscillations in the radial di- 
rection. By comparing Fig.[3ja) with Fig.[2jh), one finds 



In Fig. [4j we report the three components of the 
spin vector at Aso = F g{N — 1) = AOhcu^/a^ and 
gtl/g = 1-1- The transverse spin texture is shown in 
Fig. 5a by arrows, with color and length representing 
the orientation and the magnitude of the transverse spin 
vector (S x ^Sy), respectively. It is readily seen that the 
spin vector spirals in space and form a skyrmion-type 
texture. Quantitatively, this is most clearly illustrated in 
Fig. 5b, where we plot the skyrmion density. 
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C. Solutions of Bogoliubov equations 



Given the wave-function of the half-quantum vor- 



tex state, [(/>-(■ (p) , <j>± (p) e %v \ /y/2Tr, we now turn to 
consider its collective excitations, as described by 
the coupled Bogoliubov equations As a re- 

sult of rotational symmetry, it is easy to see 
that, the Bogoliubov wave-functions have a good az- 
imuthal quantum number to and can be written as, 
[t*t (p) , u± (p) e**, v t (p) , v± ( P ) e-^fe^/V^- There- 
fore, we have 
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To solve the Bogoliubov equation, as before we expand 
the wave-functions using 2D harmonic oscillator basis, 



u T (p) 
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This leads to a secular matrix of "Heog, whose ele- 
ments can be calculated directly using the 2D har- 
monic oscillator basis. We note that, to obtain the Bo- 
goliubov quasiparticles we cannot diagonalize directly 
the secular matrix, because of the minus sign before 
V-f (p) and v± (p) at the right-hand side of Eq. ( 23 1 . 



Instead, we should diagonalize a non-symmetric ma- 
trix Diag{+1, +1, — 1, — lj^Bog an d normalize the quasi- 
particle wave- functions according to pdp[u 2 + u 2 — 
u| — u|] = 1- The number of resulting eigenvalues is 
two times the number that we want. There are two 
branches of eigenvalues, one is positive and the other 
negative, as a result of the duality between the so- 
lution [it t (r) , u x (r) , w T (r) , v i (r)] T (with energy +hu)) 
and [vl (r) , vt (r) , u| (r) , uf (r)] T (with energy — fuS). 



We should take the positive branch. We note also that 
the Bogoliubov quasi-particles at a negative azimuthal 
quantum number m may be obtained from the negative 
branch of the solution with m > 0, because of the duality. 



1. Breathing modes 



In the case of the breathing mode (m = 0) , where 



^Boe — 



C G P+U U 

u C GP +U 



(31) 



we may have an alternative way to solve the Bogoliubov 
equation, following Hutchinson, Zaremba, and Griffin 
(HZG) [25] . By denoting collectively u = [u^ (p) , uj, (p)] 
and v = [v-f (p) , v^ (p)], we have, 

(£ G p + 2W) (u + v) = huj (u-v), (32) 
£gp (u — v) = huj (u + v) . (33) 

Let us now expand the wave-functions u ± v in terms 
of the eigenfunctions ip a of C G p with energy e a (i.e., 



a#0 

E 



1/2 
hid 



-4>a- 



(34) 



(35) 



Here, the lowest eigenstate of £qp with zero energy 
should be removed, as it corresponds exactly to the con- 
densate mode. It is easy to see that (£ G p + 2U)£>gp( u ~ 
v) = {hu) 2 (u-v) and £ G p(£ GP + ZU)(u + v) = (huj) 2 (u+ 
v). Inserting the expansion of u — v or u + v, one finds 
the secular equation, 



Y, + 2e 1 J 2 U a0 e 1 / 2 } c p = {hu) 2 



where 



pdp^l {p)U^p (p). 



(36) 



(37) 



By diagonalizing the secular matrix, one obtains the 
mode frequency uj and the coefficients c a . The latter 
should be normalized as ^ Q c 2 a = hcu, in accord with the 
normalization condition for u and v. 

We have checked numerically that the HZG so- 
lution leads to exactly the same result as the di- 
rect diagonalization of the non-symmetric matrix 
Diag{+1, +1, — 1, — l}"%Bog) if we discard the zero- 
frequency condensate mode in the latter method. 



D. Collective excitations 

In Fig. 6, we report the breathing (to = 0) and the 
dipole mode (m = ±1) frequencies as a function of the 
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Figure 6: (color online). The mode frequency of breathing 
(m = 0) and dipole (m = ±1) modes as a function of in- 
teraction strength at a fixed SO coupling Xso = 1 and at 
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Figure 7: (color online). The mode frequency of breathing 
(m = 0) and dipole (m = ±1) modes as a function of SO 
coupling at a fixed interaction strength g(N— 1) = 40fo<,x/aj_ 
and at = l.lg. 



interaction strength. With increasing interaction, the 
mode frequency decreases and seems to saturate at suffi- 
ciently large interactions. This may be anticipated from 
the point of view of two-fluid hydrodynamic behavior in 
the Thomas-Fermi regime. In Fig. 7, we report the de- 
pendence of the mode frequencies on SO coupling. In 
the absence of SO coupling, the breathing mode with 
lj = 2ujj_ and the dipole mode with u = cjj_ are the 
exact solutions of quantum many-body systems in har- 
monic traps. At a finite SO coupling, however, we find 
that these two solutions are no longer exact. The rela- 
tive deviations of the breathing mode and dipole mode 
a t ^so = 1 are about 10% and 30%, respectively. 

In Fig. 8, we plot the Bogoliubov wave- functions of 
the lowest four breathing modes at Xso = 1- g{N — 1) = 
AOhuj^/a^ and — l.lg. We find that the density re- 
sponse is mainly carried by u^(p) and u±(p) components. 
With increasing mode frequency, more and more nodes 
appear in u^(p) and u^(p). In contrast, the response in 
v-i-(p) and Vi(p) is relatively weak and the curve shape is 
nearly unchanged as the mode frequency increases. 



Figure 8: (color online). Bogoliubov wave- functions of the 
lowest four breathing modes at Xso = 1, g(N — 1) = 
40huj±/a 2 L and gf± = l.lg. The mode frequencies are indi- 
cated in Fig. 7b by solid symbols. 



E. Dynamical Calculations 

To investigate the dynamical properties of the sys- 
tem, we also perform direct simulations of the system 
by real-time propagation of the ground state under per- 
turbation. To do this, firstly we obtain the ground state 
by solving the coupled GP equations in Eqn. ([3| using 
the TSSP technique. The half-quantum vortex ground 
state is perturbed in various ways. We observe that the 
mode frequencies obtained by dynamical simulation agree 
well with those obtained by solving Bogoliubov equations 
(shown in Fig. [6]). 

Breathing mode analysis, m = 0: We excite the 
monopole mode by weak relaxation of the trapping fre- 
quency at time t = 0, and letting the system propagate in 
real-time. As the breathing mode excitation is isotropic 
in x-y space, it is sufficient to observe the dynamic re- 
sponse of the collective coordinate along one axis, say, 
the x-axis. Here, we pick the mean square of the center- 
of-mass coordinate as the quantity of interest: 

2 _ J \(j) a \ 2 x 2 dxdy 
{X)a fl^dxdy ' 

where a = f, 4,-spin components. In Fig. [9] (a),(b), we 
plot the time response of (x (t)) a for a typical param- 
eter set. In Fig. [9] (c),(d), we show the correspond- 
ing frequency response by plotting the single-sided am- 
plitude spectrum |(x 2 (u;))| (T , which are just the Fourier 
transforms of (x 2 (t)) a . We observe frequency peaks at 
~ 0.46, 1.8, 2.18 and at 3.40 (not shown). We note 
that these values exactly match with the mode frequen- 
cies obtained for this parameter set by solving Bogoliubov 
equations, shown in Fig. [6^b). 

Dynamical calculations also reveal the coupling be- 
tween the center-of-mass motion and the internal spin de- 
grees of freedom, a trademark signature of spin-orbit cou- 
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Figure 9: (color online). (a),(b): Dynamic response of the 
mean square of the center-of-mass coordinate in s-direction of 
f- and J,- spin components respectively. We have shifted the 
curves by subtracting the time-averaged {x 2 {t)) a . Without 
this shift, the Fourier spectrum as shown in (c) and (d) will be 
dominated by a large peak at cj = 0. (c),(d): Corresponding 
single-sided amplitude spectrum of the collective coordinate. 
Parameters used: Xso = 1-0, g(N — 1) = AOhuj^/a'^, gf±/g = 
1.1. 



pled systems. We shall now discuss the dynamic response 
of the population difference An = J dr (\<fi^\ 2 — |^| 2 )- In 
Fig. 10 'a), we plot the time response of An(i) for the 
same parameter set mentioned in Fig. [9] In Fig. [lO^b), 
we show the corresponding frequency response by plot- 
ting the single-sided amplitude spectrum |An(w)|. We 
observe frequency peaks at ui/ui± ~ 0.46, 1.8, 2.18 and at 
3.40 (not shown), exactly matching with the mode fre- 
quencies obtained in Fig. [9] This analysis clearly shows 
that the population transfer between the two spin com- 
ponents shares a similar dynamic response with the col- 
lective motional coordinate. In this aspect, response of 
An in a spin-orbit coupled spinor BEC (shown here) is 
similar to the effects observed in the presence of internal 
Josephson coupling in multi-component condensates |27| . 
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Figure 10: (color online), (a) Dynamic response and (b) 
single-sided amplitude spectrum of population difference An 
for the same parameter set used in Fig. [9] 

Dipole mode analysis, m = ±1: We excite the dipole 
modes by displacing the trap in x-direction by a small 
amount at time t = 0, and letting the system propagate 
in real-time. We observe the dynamic response of the 



{x) a 



j \(t><j\ 2 xdxdy 
J \<f> a \ 2 dxdy 



plot 



the time response of 
x-direction of 'f- and |- 



In Fig. [TT] (a),(b), we 
this collective coordinate in 
spin components for a typical parameter set. In 
Fig. 11 (c),(d), we show the corresponding frequency 



response by plotting the single-sided amplitude spec- 
trum Kx^))!^. We observe frequency peaks at uj/ujj_ ~ 
0.05,0.43,0.70,1.25,1.34, (shown) and at 2.5, 2.64, 2.76 
(not shown). We note that these values exactly agree 
the mode frequencies obtained for this parameter set by 
solving Bogoliubov equations, shown in Fig. [6^ a), (c). 



In the inset of Fig. 11 a), we show the dynamics of 
the center-of-mass coordinate. It is important to note 
that even though the trap is displaced only in the x- 
direction, we also observe a similar dynamic response in 
y-direction of both spin components (only f-spin compo- 
nent shown). This behavior occurs due to the vorticity 
induced by the spin-orbit coupling — the vortex state 
experiences a Magnus force that is perpendicular to its 
motion. Hence a displacement in the x-direction induces 
a motion along the y-direction. Furthermore, the trace of 
the center-of-mass and its magnitude are affected by the 
strength of the inter-particle interactions and the spin- 
orbit coupling induced population transfer, as observed 
in the case of the breathing mode excitation, between the 
f- and J.- spin components. 
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Figure 11: (color online). Parameters used: Xso = 1.0, g(N — 
1) = 40huj± /a\ , i = 1.1 g. (a),(b): Dynamic response of 
the center-of-mass coordinate in x-direction of f- and J,- spin 
components respectively. The inset in (a) shows the dynamics 
of the center-of-mass coordinate over 12 trap periods and the 
filled (red) marker denotes the initial position. (c),(d): Cor- 
responding single-sided amplitude spectrum of the collective 
coordinate. 
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IV. INSTABILITY ANALYSIS AND PHASE 
DIAGRAM 

We are now ready to analyze the parameter space for 
the existence of half-quantum vortex state. It will be- 
come unstable with respect to increasing the interaction 
strength or decreasing the ratio g^i/g. The instability 
could be indicated from some energy considerations and 
from the softening of collective density modes. 



A. Superposition instability 



As we mentioned earlier, for any half-quantum vor- 
tex state, </>(r) = [(f)^(p), (f>i(p)e lv ] T /V^n, we always 
have a degenerate time-reversal partner state, T0(r) = 
[4>\ y {p)e~ lip , — c/)^(p)] T /V2n. There is an instability for 
half-quantum vortex state with respect to a superposi- 
tion state, which with equal weight takes the form, 



Mr) 



1 

V47T 



MP) + Mp)e- l(v - Vo) 



(38) 



Here (fo is an arbitrary azimuthal angle. The energy 
difference between the superposition state and the half- 
quantum vortex state is given by, 



ae gp = (^-y-i) ^ (r)] . 



(39) 



Therefore, if W^[(/)(r)] > 0, the half-quantum vortex state 
is stable only when g < g^±. 
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Figure 12: (color online), (a) The W^-function as a function of 
SO coupling at g(N — 1) = huj±/a 2 L and g^i = l.lg. (b) The 
IV-function as a function of interaction strength at Xso = 
1 and g-\\, = l.lg. (c) The instability of the lowest dipole 
mode frequency cj m= _i with decreasing g^i/g at Xso = 1 
and g{N - 1) = 20hco ± /a 2 L . 



In Figs. 12 a) and (b), we check the VF-function of the 



interactions. It always appears to be positive, though 
the interactions tend to decrease its absolute magnitude. 
Hence, there must be a quantum phase transition occur- 
ring at the isotropic point g = g-j-j,. Once g > g^, a 
superposition state with density pattern, 



"t.4- 



1 

27 



± 0f^4. cos (ip - ip ) 



(40) 



becomes preferable. The 2D contour plot of this density 
pattern with tpo — is schematically shown in the inset 
of Fig. 1 (in the phase IIA). 

In general, in passing the quantum phase transition 
point, we would observe softening of a particular mode 
frequency. As the superposition state involves a time- 
reversal state with angular momentum m = —1, the 
lowest dipole mode with m = — 1 may become unsta- 
ble. In Fig. [l2"{c), we plot the lowest dipole mode fre- 
quency w m= _i as a function of g-\^/g at Xso = 1 & n d 
g(N — 1) = 2§huj_\_/d]_. Indeed, with decreasing g^i/g, 
the mode frequency w m= _i decreases and approaches to 
zero exactly at the phase transition point. 



B. Instability to high-order angular momentum 
components 

There is another instability for the half-quantum vor- 
tex state, occurring with increasing the interatomic in- 
teractions. With sufficiently large interactions, we an- 
ticipate that the state with high-order azimuthal angu- 
lar momentum will energetically become favorable. For 
example, let us consider a condensate state with an az- 
imuthal angular momentum m = 1 (the 3 2-quantum 
vortex state), which has the form, 



*=i (r) 



1 



Mpy v 
Mpy 2v 



(41) 



The GP energy of this state can be obtained by solv- 
ing the GPE equation as before, except that we need 
to take Rki (p) and Rk2 (p) as the expansion functions 
for $t(/?) and 4>\{p), respectively. Its degenerate time- 
reversal partner state has an azimuthal angular momen- 
tum m = —2. 



half-quantum vortex state in the presence of interatomic 



It is easy to see from Fig. 13 a) that beyond a critical 
interaction strength the condensate state with m = 1, 
</> m= i(r), is lower in energy than the half-quantum vor- 
tex state, cf> m= o(r). We note, however, that the critical 
interaction strength determined in this way is not accu- 
rate, as a superposition state of </> m= o(r) and </> m= i(r) 
may already become energetically more preferable than 
4> m =i(j) at a smaller interaction strength. 

An accurate determination of the threshold could be 
obtained by monitoring the instability in a particular 
collective mode. As the condensate state may preserve 
a well-defined parity, we find that the instability oc- 
curs in the lowest quadrupole mode with m = — 2. In 
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Figure 13: (color online), (a) GP energy of the 3/2-quantum 
vortex state m=1 (r) and of the half-quantum vortex state 
m=o (r) as a function of interaction strength at Xso = 2 
and g-fi/g = 1.1. Beyond a critical interaction strength as 
indicated by an arrow, m=1 (r) becomes energetically favor- 
able, (b) The corresponding lowest quadrupole mode fre- 
quency u) m =-2- It becomes unstable beyond a threshold g c . 




Figure 14: (color online). Phase diagram at g^i = g and 
g-fl — 2g. The critical interaction strength has been shown as 
a function of SO coupling. 



Fig. 13 'b), we report the lowest quadrupole mode fre- 
quency ui m= -2 a s a function of the interaction strength. 
As the interaction increases, the real part of mode fre- 
quency decreases down to zero and then, the imaginary 
part becomes positive, indicating clearly that this mode 
will exponentially grow if the condensate is initially in the 
half-quantum vortex configuration. The condensate then 



starts to involve high-order angular momentum compo- 
nents. The critical interaction strength g c can be simply 
determined from the softening of the mode frequency, 
u m =-2(g = 9c) = 0. 

In Fig. [14] we present critical interacting strength as 
a function of SO coupling at = g and g^ = 2g. The 
solid line at the isotropic point g^^/g has been recently 
calculated by Xiang-Fa Zhou and Congjun Wu by using 
an imaginary time evolution method [P31I2B] . Our results 
are in excellent agreement with theirs. We find that at 
smaller SO coupling the critical interaction strength de- 
creases rapidly with increasing g^/ ' g. 



C. Instability against anisotropy in SO coupling 
strength 

So far we have focused our attention on the half- 
quantum vortex state supported by an isotropic 2D har- 
monic trap subject to an isotropic Rashba SO coupling. 
Here we discuss the effect of the anisotropy in SO cou- 
pling strength on the stability of half-quantum vortex 
state. The effect of the trap anisotropy will be discussed 
in the next subsection. In the context of ultracold gases, 
anisotropic Rashba spin-orbit coupled was first discussed 
in Ref. [6] and the coupled GP equations were solved for 
a many-body system in the absence of the trap and in 
the restricted scenario when g^ = g. Here, we move 
beyond these restrictions and discuss the ground state 
of the system. We write the SO coupling term in the 
form Vso = —i(X y a x d y — X x a y d x ), where X x ,X y are SO 
coupling strengths in the two perpendicular directions. 
By including this SO coupling term and solving the cou- 
pled GP equations under the Hamiltonian as given in 
Eq. (J3j) using the TSSP technique, we obtain the ground 
state wavefunction at various values of anisotropy in SO 
coupling represented by X x /X y . 



In Fig. 15 we plot the 
corresponding ground state density profiles of 4,-spin com- 
ponent for an SO coupling strength of X x = 4.0, and for 
various values of X x /X y . 

We see from Fig. [TB^a) that the half-quantum vortex 
state is indeed the ground state (already mentioned in 
Fig.JIJb)) for the parameter set: g(N — 1) = 0.1hu)±/a\, 
gfl/g — 1.1, X x = 4.0 and X x /X y — 1.0. We shall now 
analyze the pattern in which the density profile changes 
with anisotropy in SO coupling strength as shown in 
Fig. 15^b)-(d). It is evident from the density distribu- 
tions in Fig. |15| that the half-quantum vortex state is 
unstable even against small anisotropy in SO coupling 
strength. Adopting a similar method as presented in 
Ref. |29| . we analyze this systematically by expanding 
the wavefunction of ^--component in an orthogonal ba- 
sis set of the form: = £„ f n (p) (2ll+1 ) v , where 
n measures the vorticity, and f n {p) absorbs the nth 
mode's contribution in radial direction. We quantify the 
weights of the wavefunction in t he n th mode by comput- 
ing a n = J dp\f n (p)\ 2 . In Fig. 16 we plot the weights 



a n relative to ap computed for half-quantum vortex state 
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D. Instability to anisotropy in trap potential 




Now we examine the effect of anisotropy in the trap- 
ping potential, but with isotropic SO coupling, on the 
stability of half-quantum vortex state. We write the trap- 
ping potential in the form V(x, y) = M(u> 2 x 2 +uj 2 y 2 )/2 = 
MlJ\{x 2 + f 2 y 2 )/2, where uj x = w±,u) y — f y oj± are trap- 
ping frequencies in x- and y-directions respectively. We 
again obtain the ground state wavefunctions at various 
values of f y by solving the coupled GP equations using 
the TSSP technique. In Fig. [IT] we plot the correspond- 
ing ground state density profiles of |-spin component for 
an SO coupling strength of Xso — 4.0, and for various 
values of trap anisotropy ranging from to 10%. 






Figure 15: (Color online) Plot of the ground state density 
profiles of ^-spin component for the parameter set: g(N—l) = 
0.1hu)±/a 2 Ll g-fi/g = 1.1, \x = 4.0, but with varying ratios of 
\ x /\ y . (a) Isotropic case: \ x /\ v = 1.0, (b) X x /X y = 1.01, 
(c) X x /\y — 1.05, (d) Xx/Xy — 1.1. Viewing angle is slightly 
tilted for aesthetic purposes. 



y^o 



with X x /X y — 1.0. As we would expect, for this isotropic 
case, do = 1 and a n — for n / 0. As anisotropy in 
SO coupling strength increases, more and more n ^ 
components will be mixed into the ground state. 
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Figure 16: (Color online) Plot of the weights of ground-state 
wavefunction of ^-spin component - corresponding to the den- 
sity profiles in Fig.[l5]- in the nth mode. The weights are nor- 
malized with respect to ao computed for half-quantum vortex 
state with X x /X y — 1.0. (a) Isotropic case: \ X /X y = 1.0, (b) 
X x /X y = 1.01, (c) X x /X v = 1.05, (d) Xx/Xy = 1.1. 
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Figure 17: (Color online) Plot of the ground state density 
profiles of ^-spin component for the parameter set: Xso = 4.0, 
g(N — 1) = 0.1ftajj_/a^, g^i/g = 1.1, but with varying ratios 
of f y — Lo y /uj x . (a) Isotropic case: f y — 1.0, (b) f y — 1.01, 
( c ) fy = 1-05, (d) f y — 1.1. Viewing angle is slightly tilted 
for aesthetic purposes. 

We see from Fig. [l7[a) that the half-quantum vortex 
state is indeed the ground state (already mentioned in 
Fig. 15 a)) for the parameter set: Xso — 4.0, g(N — 1) = 
O.lhu^/a 2 ^, g\i/g = 1.1. We shall now analyze the 
pattern in which the density profile changes with trap 
anisotropy Fig. |17[ b)-(d). It is evident from the density 
distributions in Fig. |17| that the vortex core becomes 
increasingly anisotropic with increasing f y . We analyze 
this systematically by expanding the wavefunction of J,- 
component in an orthogonal basis set and quantifying the 
weights in the nth mode by a n , as mentioned in Sec. |IV C| 
In Fig. 18 we plot the weights a n relative to ao com- 



puted for half-quantum vortex state with /„ = 1.0. As 



we would expect, 
ao = 1 and a n = 



for the isotropic case with f y = 1.0, 
for n 7^ 0. As trap anisotropy in- 
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creases, we observe that the ground state is a mixture 
of n 7^ components as well. Nevertheless, we see that 
the trap anisotropy has a much smaller effect on the half- 
quantum vortex state than the anisotropy in the SO cou- 
pling strength. 
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Figure 18: (Color online) Plot of the weights of ground-state 
wavefunction of ^-spin component - corresponding to the den- 
sity profiles in Fig. [17]- in the nth mode. The weights are nor- 
malized with respect to ao computed for half-quantum vor- 
tex state with f y = 1.0. (a) Isotropic case: f y — 1.0, (b) 
/„ = 1.01, (c) /„ = 1.05, (d) /„ = 1.1. 



We have found that: 

(1) The condensate is in a half-quantum vortex state, 
if the intra-species interaction g is smaller than inter- 
species interaction and, if the interaction strength is 
below a threshold g c . We have calculated the threshold 
by monitoring the unstable quadrupole mode with an 
azimuthal angular momentum m — — 2. A phase diagram 
for the half-quantum vortex state is therefore determined, 
as given in Figs. fl] and [L4| 

(2) The half-quantum vortex state (the phase I) will 
turn into a superposition of two degenerate half-quantum 
vortex states (the phase IIA) if g > and will start to 
involve high-order angular momentum components (the 
phase IIB) if g > g Cl where g c depends critically on the 
ratio g-\\./g. The half-quantum vortex state is unstable 
against small anisotropy in SO coupling strength and 
large anisotropy in trapping potential. The state tends to 
be a superposition of higher angular momentum states. 

(3) In the presence of spin-orbit coupling, the behavior 
of collective density modes becomes complicated. In par- 
ticular, the breathing mode with uj = 2uj± and the dipole 
mode with uj = lj±_ are no longer the exact solutions of 
the many-body system. 

(4) The condensate wave-functions in the phases IIA 
and IIB are yet to be determined using the time-splitting 
spectral method for GPE. These wave-functions break 
the rotational symmetry. We anticipate that interesting 
density patterns will emerge in the limit of very large 
interatomic interactions. This is to be explored in future 
studies. 



V. CONCLUSIONS 

In summary, we have investigated systematically the 
ground condensate state of a spin-orbit coupled spin- 
1/2 Bose gas confined in two-dimensional harmonic 
traps. The density distributions and collective density 
excitations have been obtained respectively by solving 
the Gross-Pitaevskii equation and Bogoliubov equation, 
which are generalized to include the spin-orbit coupling. 
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